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Abstract 

Ab initio electronic structure methods give accurate results for small systems, but do not scale 
well to large systems. Chemical insight tells us that molecular functional groups will behave ap- 
proximately the same way in all molecules, large or small. This molecular similarity is exploited 
in semiempirical methods, which couple simple electronic structure theories with parameters for 
the transferable characteristics of functional groups. We propose that high-level calculations on 
small molecules provide a rich source of parametrization data. In principle, we can select a func- 
tional group, generate a large amount of ab initio data on the group in various small-molecule 
environments, and "mine" this data to build a sophisticated model for the group's behavior in 
large environments. This work details such a model for electron correlation: a semiempirical, 
subsystem-based correlation functional that predicts a subsystem's two-electron density matrix as 
a functional of its one-electron density matrix. This model is demonstrated on two small systems: 
chains of linear, minimal-basis (H-H)5, treated as a sum of four overlapping (H-H)2 subsystems; 
and the aldehyde group of a set of HOC-R molecules. The results provide an initial demonstration 
of the feasibility of the approach. 
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I. INTRODUCTION 



Canonical ab initio electronic structure methods provide highly accurate electronic struc- 
tures for small systems of (9(10) atoms. However, these methods are too computationally 
intensive to apply to large systems. The formal scaling of computational effort for ab initio 
calculations on an A^-electron system ranges from 0{N^) for Hartree theory, to 0{N^) for 



MP2, to (9(e^) for the exact, full-configuration-interaction (fuU-CI) solution Ab initio 
("first principles") calculations always begin with a minimal amount of information about 
the system (e.g. an initial geometry and a basis set), determining practically all of the 
system's features at runtime. 

The computational effort of ab initio calculations can be mitigated using two physically- 
motivated approximations: 0(N) and semiempirical approximations. 

0{N) approximations are based on the principle of nearsightedness Q, which states that 
the interactions between parts of a molecule are largely local in character. (A discussion of 



nearsightedness can be found in Ref. 



0(N) approximations have been developed for 



every part of an ab initio calculation, from fast multipole methods for Coulomb effects 



to divide-and-conquer ^ and other 




d (SCF) 
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calculations, to treatments of electron correlation 
schematic of a nearsightedness-based approximation as outlined in Ref. j2| is shown in Fig. 

Semiempirical approximations are based on the principle of molecular similarity: that 
the properties of atoms and functional groups are largely conserved in different molecules. 
This principle formalizes the chemical insights that methyl groups are relatively small and 
nonpolar, halides are electron- withdrawing, and so forth. Ab initio calculations spend much 
of their time in re-calculating the transferable characteristics of functional groups. Semiem- 
pirical approximations replace the ab initio Hamiltonian with a simpler model Hamiltonian, 
which contains parameters that capture the transferable characteristics of functional groups. 
Examples of these parameters include force constants in molecular mechanics 
tonian matrix elements in semiempirical quantum-mechanical approximations [2 

One of the benefits of 0{N) methods is their controllability. 0(N) approximations yield 
well-defined changes in accuracy and computational effort. The decision to use an (9(N) 
approximation can be made a priori based on the size of the system of interest P]. 

Unfortunately, semiempirical methods usually involve a significant trade-off between com- 



1911 or Hamil- 
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putational effort and accuracy. Semiempirical methods are much less accurate than ab ini- 
tio methods for many systems. This has led to the widespread use of "hybrid" QM/MM 
methods, a nearsightedness-based tradeoff between ab initio accuracy and semiempirical 
speed (22, |23| . Our goal is to systematically improve semiempirical theory. 

Most existing semiempirical methods are based on models that were designed to 
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semiempirical methods are now 



27|), we believe that the existing 



be parametrized to experimental data. Thou^ 
parametrized using ab initio results (e.g Refs. \24 . 
methods may not take full advantage of the possibilities inherent in ab initio parametrization. 
Ab initio calculations on small molecules can give orders of magnitude more parametriza- 
tion data than can be readily obtained from experiment. They also yield information that 
is more directly relevant to a semiempirical model's parameters. 

Nearsightedness and molecular similarity suggest that we can model large systems as 
the sum of contributions from different functional groups. This implies that a sufficiently 
rich data set of a functional group in small molecules will contain all information needed 
to describe the functional group in molecules of arbitrary size. Our overall approach is to 
generate rich data sets on the behavior of functional groups by doing a large number of 
highly accurate ab initio calculations on the group in a set of small-molecule environments. 
This paper investigates whether a semiempirical model parametrized to this sort of small- 
molecule data can give ab initio accuracy for larger molecules. 

This approach is fairly general. It requires only that the semiempirical model can describe 
a system as a sum of subsystem contributions. For example, a semiempirical model that 
predicts the amplitudes of delocalized wavefunctions would not be compatible with this 
approach. 

The current work details our first implementation of this approach: a semiempirical 
subsystem-based treatment of electron correlation. We model the system in terms of its one- 
and two-electron density matrices in an atomic orbital basis set (Sec. Ill A|) . Subsystem two- 
electron density matrices are combined to model the two-electron density matrix of the entire 
system. This model was chosen because it treats an important problem in contemporary 
electron structure theory (electron correlation), and because the predicted outputs (two- 
electron density matrices) are much easier to obtain from ab initio calculation than from 
experiment. 
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II. METHODS 



A. Semiempirical model for electron correlation 

Our semiempirical model treats electron correlation by predicting subsystem two-electron 
density matrices as a functional of subsystem one-electron density matrices. A system's 
one- and two-electron density matrices ^D, are obtained from its normalized A^-electron 
wavef unction |$) as 

^D{a,h) = maia,\^) (1) 

2L)(ac,6rf) = 1/2 ($ I a^ajabad I $) (2) 

in second quantization with one-electron basis functions {|0a)}- For an N-electron system, 
the trace of equals N and the trace of equals the number of unique electron pairs, 
1/2 A^(A^ — 1). The electron-electron interaction energy of a system {E2) is obtained as the 
trace over the product of the two-electron integrals and the two-electron density matrix 

^2 = ^ {ac\bd) ^D{ac,bd) (3) 

abed 

(ac\hd) = [ (iri(ir20*(ri)0*(r2)-^ ^ r0fe(ri)0d(r2) 

J |r2-ri| 

The electron density in real .pace i. the diagonal ot the one-electron density matrix: >fl(r) . 

($1 ajor 1$) [28]. and provide a complete description of a system whose Hamiltonian 

ir~i 

contains only one- and two-body interactions 29]. 



The two-e 



ectron density matrix obtained from 1$) can be expressed as a cumulant 



expansion j30|, |31| 

'^D{ac,bd) = 1/2 ^D{a,b) ^D{c,d) (4) 
- 1/2 ^D{a,d) ^D{b,c) 
+ ^A{ac,bd) 

where the three terms on the right-hand side of Eq. 0] are denoted Coulomb, exchange, 
and correlation contributions to ^D. The connected pair density matrix is that part 
of that cannot be written as a simple function of ^D. The Coulomb and exchange 
contributions to in Eq. |3 are well-approximated at the Hartree and Hartree-Fock levels 
of theory, respectively. However, accurate ah initio treatment of the connected pair density 
matrix requires expensive high-level methods. 



Density functional theory (DFT) is a formally exact method for treatiiig a system of 



29 



32, |3i 



The 



interacting electrons exclusively in terms of its one-electron density 
electron-electron interaction energy of Eq. |21 is treated as the sum of a Coulomb term and 
an exchange-correlation correction Exc^ such that the electrons move in a potential that is 
corrected by the exchange-correlation potential fxc(i") = (^{Exc)/^(^D{r)) . DFT is im- 
plemented by approximating vxc as a functional of electron density: vxc = vxc[^D{r)]. 
(In Kohn-Sham DFT, the kinetic energy is decomposed into the kinetic energy of the 
Kohn-Sham orbitals plus a den sity -dependent correction, which is incorporated into vxc 
via e.g. adiabatic connection [2^. Our Vcorr functionals (Eq. ^ do not include a ki- 
netic energy correction.) Following Eq. |21 and Eq. 0J Exc may be obtained as the 
trace over the exchange and correlation contributions to the two-electron density matrix: 
Exc = J2 {(^c\bd) {—1 / 2 ^D{a,d) ^D{b,c)+ ^A(ac, foe/)). Thus, the correlation component 
Vcorri^D] of a systcm's exchange-correlation functional can be obtained from the first deriva- 
tive of a functional that predicts a system's connected pair density matrix as a function 
of its electron density matrix ^D: 



Vcorr[^D]{a' ,b') = {ac\bd) 



abed 



(5) 



Explicit treatments of vxci^D] in terms of the two-electron density include various analyses 



x5(^A['D] (ac, bd)) /5CD{a,b')) 



of the real-space exchange-correlation hole 



35 



In this work, we define the correlation energy Ecorr as the expectation value of the con- 
nected pair density matrix: Ecorr = I] (ac|6(i) ^A(ac, Correlation energy can also be 
defined as the difference in energies predicted by configuration-interaction and Hartree-Fock 
calculations: Ecorr = Eci — Ehf- The latter definition includes the effects of relaxation, 
e.g. the expectation value of ^Dci — ^D^f- In contrast, the former definition yields the 
correlation energy corresponding to a single choice of ^D, and is therefore consistent with 
the definition of Ecorr used in DFT and MP2 calculations. 

Both and ^A can be treated using the nearsightedness approximation. Several 
nearsightedness-based treatments of exist, includingdivide-and-conquer methods that 
partition into subsystem contributions as in Fig. ^ Nearsightedness- 
based treatments of ^A include the CXN) treatments of electron correlation cited previ- 
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We recently developed the "localized reduced 



density matrix" (LRDM) method [4l|, a divide-and-conquer style treatment of ^A. LRDM 
assembles a large system's atomic-orbital-basis ^ A from the results of ab initio calculations on 
overlapping subsystems. Like other divide-and-conquer methods, LRDM is non-variational. 

In the current work, we use LRDM as a framework for semiempirical subsystem-based 
approximations for DFT correlation functionals Vcorr[^D]. We generate semiempirical func- 
tionals that predict the matrix elements of a subsystem's ^A as a function of the subsystem 
electron density matrix ^D: ^A [^D] (see Eq. [7| and Eq. IHl below). An approximate ^A 
for a large system is obtained by combining subsystem ^A [^D] predictions using LRDM. 
Our results indicate that basis-set ^A [^D] functionals can provide good results for multiple 
subsystem geometries fSec. Illl|) . 

LRDM can treat long-range correlations (dispersion interactions) by d oing ab initio cal- 



4l|. In the cur- 



culations that include correlation in two disjoint regions of a molecule 
rent work, we do not model these long-range interactions. Therefore, our subsystem-based 
Vcorri^D] functionals, like standard DFT Vcorr functionals 0, cannot treat dispersion 
interactions. 



Our subsystem-based Vcorr[^D] functionals are very different than the standard DFT 

331 13, 3|- Other groups have 

developed vxcl'^D] functionals that are semiempirical 



functionals derived from the homogeneous electron gas 



fitted to high-accuracy ab initio data 
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47[ , subsystem-based j48|, |49| or 



but to our knowledge the current method is 
unique in combining semiempirical methods with a nearsighted, molecular-similarity-based 
treatment of ^A. 



B. Parametrization method 

Our approach is to develop semiempirical models that are parametrized using rich data 
sets of small-molecule ab initio calculations. These rich data sets allow us to use data mining 
methods in the parametrization stage. "Data mining" refers to computational methods for 
analyzing large data sets and automatically extracting previously unknown dependencies 
between the data js^ . Other data-mining treatments of electron correlation include a neural- 
network exchange-correlation potential fitted to data from many molecules 0| , and a model 



for the correlation energy between pairs of widely separated, localized electrons 54 1. 



Data mining methods can determine two types of relationships between data. The first 
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is the system's dimensionality: which input and output descriptors are most important for 
describing the data set. (Here, "descriptor" is a generic term for a type of data used by a 
model. For example, the input and output descriptors of our [^D] functionals are ^D{a, h) 
and ^A(ac, hd) matrix elements.) The second type of relationship that can be determined 
by data mining is the functional form of the [input descriptor] —> [output descriptor] relation. 
In the current work, we assume a quadratic input-output relation and focus on dimensional 
reduction. 

A flowchart of the data-mining process for a functional group is as follows. 

1. Choose an initial set of input and output descriptors, and a fit function to relate 
them. As discussed in Sec. HI subsystem-based models require input and output de- 
scriptors that describe electronic structure in terms of local information (e.g. electron 
densities). Since the models are meant for use within semiempirical models, the in- 
put descriptors should be obtainable from a simple approximate Hamiltonian (e.g. the 
DFT Hamiltonian). The fit function can be anything from a polynomial fit to a neural 
network. 

2. Generate an initial data set of ah initio calculations on the functional group in var- 
ious small molecules and environments. Extract the values of all input and output 
descriptors for each point in the data set. 

3. Split the data set into training and testing subsets. 

4. Reduce the dimensionality of the data set, by using (for example) principal component 
analysis to determine a few combinations of descriptors that capture most of the 
variation in the data set. The model will be parametrized on this dimensionally- 
reduced input and output data. 

5. Parametrize the model using the training subset of the small-molecule data. 

6. Test the model on the testing subset of the small- molecule data, and on larger 
molecules. 

As stated above, our initial focus is on the dimensional reduction of and (step^. 
For a system with M basis functions, and ^A have 1/2 (M^ + M) - 1 and 1/8 (M^ + 
2M^ + 5M^ -|- 4M)-1 degrees of freedom, respectively. Without dimensional reduction, even 
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a small subsystem (e.g. M = C(10)) has far too many output degrees of freedom for a 
[^D] functional to be useful. We use principal component analysis (PCA) to decompose 
and into a set of importance-weighted basis functions. To illustrate, PCA on a set 

of subsystem electron density matrices {^D„} describes each density matrix as 

'D,{a, b) = 'Da.g{a, b)+J2 o.i ^D,{a, b) (6) 

i 

where ^Davg is the average electron density matrix, ^Di are the principal components, and 
the standard deviation of the expansion coefficients Co-j, evaluated across the data points a, 
decreases with increasing i jssl. 

In this work, we use a quadratic function to predict the first few (most important) ^A 
components from the first few components. A quadratic function is the lowest-order 
polynomial of ^A[^D] for which the associated correlation-energy functional Vcorr[^D] (Eq. 
is not a constant. The ^A [^D] functionals fit the first C2 principal components of ^A as a 
function of the first Ci components of such that 

C2 

^A[^D]{ac, bd) = ^Kvgiac, bd) + ^{^A^^ac, bd) (7) 

i 

i 

where ^Aj are the principal components of ^A, (^D\^Di) is the projection of the argument 
one-electron density matrix onto the iih principal component 

[^D\'D,)=Y,{^D{a,b)- 'Da,g{a,b)) 'D,{a,b) (8) 

ab 

and {ajjjijjaij} are fitted parameters. Each component of the two-electron density matrix 
is fit independently of the others. The subsystem DFT correlation energy operator Vcorr[^D] 
is obtained from ^A [^D], following Eq. as 

iv,orr['D])ia',b') = J2 {ac\bd)f^{'A,iac,bd) 

abed j 

Ci 

'D,{a',b'){-fi,+2ai,{'D\'Di))} (9) 

i 

The Vcorri^D] functional of a large system is obtained by overlaying subsystem contributions 
as in LRDM (Fig. QJ. The degree of dimensional reduction can be seen by comparing the 
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number of fitted components Ci and C2 to the total number of degrees of freedom in 
and 2 A. 

The approach discussed here can be used to construct many different kinds of semiempir- 
ical model based on the choice of input and output descriptors. For example, we have begun 
work on a semiempirical model of core polarization in effective core 

input descriptors are the valence electron density matrix and one-electron Hamiltonian, and 
the output descriptors are the core electron density matrices. This work will be discussed 
in a future publication. 

The principal computational challenge of this approach is the steep scaling of the amount 
of training set data required. The [^D] functional of an M-orbital subsystem will have 
1= 0{yP) ^Di input components [56]. In general, a function with I input components 
must be parametrized using 0{e^) data points js^. Because of this, we have focused our 
initial work on proof-of-concept treatments for small model systems in minimal basis sets. 
Successful application of our method to these systems will justify expanding our focus to 
larger, chemically-interesting systems and larger basis sets. 



C. Error Decomposition 

Our ^A l^D] functionals contain three distinct approximations. The first approximation 
is nearsightedness: the connected pair density matrix ^A is assumed to be well-described 
by a decomposition into overlapping subsystems. The second approximation is that each 
subsystem ^A is assumed to be well- described by a relatively small number of principal 
components {C2 in Eq. [7j). The third approximation is that each subsystem ^A component 
is assumed to be well-described by the functional in Eq. [7| 

We can isolate the effects of each of these assumptions using three kinds of approximate 
connected pair density matrix (see Table H}. The first approximate connected pair density 
matrix is the "exact subsystem" connected pair density matrix: ^A^^s^fc. This is obtained by 
projecting the correct ^A onto the overlapping subsystems, and setting to zero all matrix 
elements that are not contained within a subsystem. The second is the principal component 
reduction connected pair density matrix: ^Ap^A- This is obtained by projecting the sub- 
system blocks of "^Axsub onto the ^A components that are fitted by the subsystem ^A [^D] 
functionals jsT^. The third is the connected pair density matrix obtained using the subsys- 
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tern [^D] functional and the correct one-electron density matrix: '^A\^Dexact\- Table H] 
summarizes the approximations used in generating these density matrices. 



III. RESULTS 



The remainder of this paper details demonstrations of our semiempirical subsystem-based 
\^D] functionals. We begin by demonstrating ^A [^D] functionals for a linear dimerized 
chain of minimal-basis hydrogen atoms (H-H)5, since the functional predictions can be read- 
ily compared to full-CI. Then, we demonstrate that a ^A \^D] functional for the aldehyde 
group, parametrized to data from a set of small HOC-R molecules, can extrapolate to R 
groups outside of the training set. All ah initio calculations were performed using a modified 
version of the GAMES S electronic structure program 



A. (H-H)2 and (H-H)5 systems 

The first system is linear minimal-basis (H-H)5. This system is treated as a sum of four 
overlapping (H-H)2 subsystems. We model its correlation energy by parametrizing a (H-H)2 
^A \^D] functional to data on isolated (H-H)2 molecules, and combining the (H-H)2 ^A ^D] 
predictions using LRDM js^. The functionals are parametrized to, and tested on, full-CI 
ah initio calculations. 

We generated data for both variable- and fixed-geometry molecules, yielding the four 
data sets in Table IHl Each molecule was electrostatically perturbed by randomly placing 
fractional charges (| charge | < 1) into a QA x QA x (molecule length -|- AA) box around the 
molecule, with a minimum point charge - atom separation of 1.2 A. Variable geometry 
systems had each bond length set randomly within the ranges in Table |n] 

The (H-H)2 ^A \^D\ functionals were parametrized using half of the (H-H)2 data as a 
training set (see item|21in the flowchart of Sec. Ill B|) . Separate functionals were parametrized 
for the variable- and fixed-geometry systems. The numbers of principal components included 
in the ^A ^D] functionals (Ci and C2 in Table ITT|l were selected to give good results for both 
'^/\\^D exact] and /\\^ D £,pt\ (see below). The principal component analyses were a significant 
dimensional reduction, as the and ^A of (H-H)2 contain 9 and 59 degrees of freedom, 
respectively. 
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B. Modeling (H-H)2 using (H-H)2 [^D] functionals 



The first test of the (H-H)2 [^D] functionals is how well they can predict the (H-H)2 
^A given the correct full-CI electron density matrix ^Df.xact- Fig. |21 plots predicted vs. real 
Ecorr for the (H-H)2 systems. Table ITTTl presents \6Ecorr\ errors averaged over the training- 
and testing-set data. Here, |(5i?corr| is the absolute error in the predicted correlation energy 
Ecorr = J2abcd {cLc\hd) '^A{ac,bd) (see Eq. El). The ^A [^D] Ecorr predictions are compared to 
MP2. 

The results are quite encouraging. The ^A [^D] functionals are better than MP2 at 
predicting the average value of the correlation energy: the mean absolute errors in Ecorr 
from MP2 are 40 and 150 times as large as the error for the ^A [^D] functionals (variable 
and fixed geometry, respectively). ^A [^D] functionals are also better than MP2 at predicting 
the variation of the correlation energy across the data set. This can be seen in Fig. |21 the 
slope of the predicted vs. real Ecorr values is very small for the MP2 predictions but close 
to 1 for the ^A [^D] functionals. Despite its low scatter, MP2 does not capture either the 
value or the variation in the correlation energy. 

The scatter in the ^A [^D] predictions for the fixed-geometry system can be reduced by 
parametrizing a ^A [^D] functional that uses more principal components. We parametrized a 
^A [^D] functional for fixed-geometry (H-H)2 that includes seven and eight ^A principal 
components (Ci = 7, C2 = 8 in Eq. [7j). This functional gives an between real and 
predicted Ecorr of 0.990, comparable to the 0.991 value for MP2 and better than the 0.890 
value in Fig. |21 For this functional, the average (standard deviation) "^Al^Dcxact] \SEcorr\ 
values are 0.11 (0.29) mH for the testing-set data. 

A comparison of the "^ApcA and "^Al^ Dcxact] errors in Table ITTTl shows that most of the 
error in the variable-geometry system is due to dimensional reduction of ^A, as the "^ApcA 
errors are almost as large as the corresponding ^A[^Dcxact] errors. In contrast, the error in 
the fixed-geometry system is more evenly partitioned between dimensional reduction of ^A 
and prediction of ^A from ^D. 

The training- and testing-set errors are reasonably close to each other, indicating that the 
functionals are not over-fitted. We tested a second measure of the predicted ^A, the sum of 
absolute errors in the predicted ^A matrix elements. These errors were fairly well-correlated 
with the \6Ecorr\ errors presented above (data not shown). 
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The results in Table ITTTl and Fig. |2]show that the constant returned by '^A[^Davg] is a 
surprisingly good approximation for the variable-geometry systems. This is encouraging, as 
it suggests that even the most primitive [^D] functional (e.g. a constant ^A) can work 
reasonably well for multiple subsystem geometries. Our ^A [^D] functionals all improve upon 
this primitive functional, as all ^A[^Dea:arf] errors are lower than the corresponding "^Al^Davg] 
errors. As expected, "^Al^Davg] predicts a constant correlation energy for the fixed-geometry 
systems (Fig. |21). 

C. Modeling (H-H)5 using (H-H)2 [^D] functionals 

The results in Fig. |21 and Table IIIII demonstrate that the (H-H)2 ^A [^D] functionals 
give good ^A predictions for (H-H)2. Given this, we investigate whether four copies of an 
(H-H)2 ^A [^D] functional, combined using LRDM, will suffice to describe correlation effects 
in (H-H)5. Using the (H-H)2 ^A [^D] functional on (H-H)5 tests whether the fundamental 
assumptions of nearsightedness and molecular similarity, and our implementation of these 
approximations, are correct for the (H-H)5 model system. Fig. Eland Table HVl present data 
for (H-H)5 systems, using the notation of Fig. |21and Table ITTTl 

The (H-H)5 results are also encouraging. Four copies of an (H-H)2 ^A [^D] functional, 
combined via LRDM, are better than MP2 at describing the mean and variation of E^orr for 
the (H-H)5 system. The mean absolute Ecorr errors for MP2 are 60 and 90 times the values 
for ^A [^D] functionals (variable and fixed geometry, respectively). Fig. El shows that our 
method does better than MP2 at capturing the variation in E^orr across the data set, with 
a predicted vs. real Ecorr whose slope is very small for MP2 but near one for our method. 

For the variable-geometry systems, the ^A [^D] functionals describe the (H-H)5 data to 
about the same level of accuracy (per atom) as the (H-H)2 data. The average (H-H)5 |(5i?corr| 
are about 10/4 = 2.5 times as large as the corresponding (H-H)2 values. For example, the 
average "^Al^D^xact] error is 1.91 mH for variable-geometry (H-H)2 and 3.43 mH for variable- 
geometry (H-H)5. 

For the fixed-geometry systems, the ^A [^D] functionals do not describe the (H-H)5 data 
to the same level of accuracy as the (H-H)2 data: the average (H-H)5 \6Ecorr\ are about five 
times as large as the corresponding (H-H)2 values. This error is not due to the subsystem 
decomposition: the average \6Ecorr\ of "^Axsub is only 0.02 mH (Table ITV]l . We suggest that 
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the long-range order in the fixed-geometry (H-H)5 leads to an intrinsic difference between the 
environments experienced by an isolated (H-H)2 vs. an (H-H)2 embedded in (H-H)5. Better 
(H-H)2 l^D] functionals for the fixed-geometry systems could perhaps be generated by 
using cyclic boundary conditions in the (H-H)2 data. Evidence for this conclusion is discussed 
in the Supporting Information. 

The predictions of a semiempirical model should not depend on the choice of training set 
used to parametrize the model. We parametrized (H-H)2 [^D] functionals using multiple 
choices of training set. Results are discussed in the Supporting Information. As expected, 
the functionals have only a weak dependence on training set choice. 

D. DFT calculations with [^D] functionals 

The above results test the ^A [^D] functional's predictions given the correct electron 
density matrix ^Dexact- However, the functionals are intended for use in density functional 
theory (Sec. Ill A|) where exact is not known in advance. We have implemented two methods 
for using the functionals. The first is DFT with exact exchange and the ^A \^D] correlation 
functional of Eq. El referred to as '^/\\^Dj:,pt\- The second method, like MP2, is a one-step 
post-Hartree-Fock prediction of Ecorr- Here, the correct electron density matrix ^D^xact is 
approximated as the Hartree-Fock electron density matrix ^Dhf, and the correlation energy 
is obtained non-self-consistently from "^AI^Dhf]- 

Table M presents \6Ecorr\ values for "^AI^Ddpt] and "^AI^Dhp] on the fixed- and variable- 
geometry (H-H)2 and (H-H)5 systems, for a single choice of training set. Predicted vs. real 
Ecorr for the (H-H)5 systems are plotted in Fig. HJ 

The "^Al^ Dpipp] calculations do a fairly good job of predicting the average and variation 
in ^A and Ecorr'- the average and standard deviations in |5ii^corr| are much better than MP2, 
and the standard deviations in I^E'corrl are generally smaller than the primitive "^Al^Davg] 
functional (see Tables ITTTl and HV |) . These results are encouraging, given that our Vcorr[^D] 
functional is a simple linear function of (see Eq. |^. The results from the (H-H)5 systems 
are especially encouraging: four identical, overlapping (H-H)2 Vcorrl^D] functionals give a 
reasonable prediction for the Vcorr[^D] of (H-H)5. Fig.lUshows that the relatively large aver- 
age errors in the '^Al^D^pp] \5Ecorr\ are mostly systematic error. Better ^A[^Ddpt] results 
could perhaps be generated using a more sophisticated (nonlinear) function (see fiowchart. 
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Sec. IIIBf) . The errors in the non-self-consistent '^A[^Dhf] calculations are somewhat higher 
than the self-consistent '^A['^Ddft] calculations. This is reasonable, especially given that 
^Dhf is not necessarily a good approximation for ^Df.xact- 

When \^D] functionals are combined with exact exchange, the ^/\\^Dhf] and 
'^/\\^D£,pt\ calculations give a fairly large, systematic under-estimate of the total energy. 
This is partly due to a difference between the obtained using fuU-CI and HF the- 
ory on minimal-basis (H-H)5. For a closed-shell N-electron system, the combined trace 
of the exchange and correlation parts of (Eq. 0]) equals —1/2 [2^. Full-CI cal- 
culations on (H-H)5 give with a trace < and an exact exchange pair density ma- 
trix '^Dx{cLC,bd) = —1/2 ^D{a,d) ^D{b,c) with a trace less than —1/2 N. Thus, all 
of the (H-H)2 ^A [^D] functionals return a ^A with a negative trace. However, ab ini- 
tio methods that return a single-determinant wavefunction (e.g. HF or KS-DFT theory) 
give a closed-shell ^Dx with a trace equal to —1/2 A^. Thus, for example, the approxi- 
mate two-electron density matrix returned by our non-self-consistent method '^D{ac,bd) = 
1/2 ^DHF{a, c) ^DHF{b, d) - 1/2 ^DHF{a, d) ^DHF{b, c) + ^A[^DHF]iac, bd) will always have 
a trace less than the correct value (1/2 A^(A^ — 1), see Sec. Ill Ajl . This leads to a system- 
atic under-estimate of the number of electron pairs in the system and the electron-electron 
interaction energy. One way to correct this is by renormalizing the exact-exchange ^Dx 
obtained from ^Dup or ^D^pT, such that the final predicted has the correct trace. This 
is analogous to the use of a fraction of exact exchange in "hybrid" DFT functionals such 
as B3LYP 0- 

This significantly improves the total energies: for example, the average 
(standard deviation) total energy error for variable-geometry (H-H)5 is 114.71 (13.50) mH 
for uncorrected Hartree-Fock calculations and -71.52 (7.26) mH and -160.88 (19.67) mH for 
'^A[^Dhf] with and without renormalization of ^Dx- 

E. Substituted aldehydes 

The assumption of molecular similarity implies that a ^A [^D] functional for the aldehyde 
group of HOC-R molecules should be able to extrapolate to R groups outside of its training 
set. We tested this assumption by parametrizing aldehyde ^A [^D] functionals using minimal- 
basis (ST0-3G) HOC-R molecules with six different R groups: H, F, OH, CH3, CI, and 
OCH3. Six different ^A [^D] functionals were generated from this data. Each was trained 
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on a data set that excluded data from one of the six R groups, and included half of the 
data from the other five groups. The functionals were tested for their ability to accurately 
model the aldehyde for both the five kinds of HOC-R molecules in the training set and the 
R group excluded from the training set. 

Details of the calculation are as follows. The ab initio data set contained 250 calculations 



for each of the six kinds of HOC-R molecules. Each calculation had random geometric [6 
and electrostatic j^| perturbations similar to those in the variable-geometry (H-H)5 chains 
above. Ab initio calculations were performed using MP2, as the different-sized HOC-R 
groups required a size-consistent method and fuU-Cl was prohibitively expensive. The alde- 
hyde [^D] functionals were fitted to the MP2 and the relaxed [62] of the aldehyde 
group. The aldehyde functionals' performance was characterized by their ability to reproduce 
the "aldehyde correlation energy" defined as E^J^^ = J2 {cic\bd) ^A(ac, bd) ; {abed} G HOC. 
All functionals used 40 and 30 ^A principal components. This was a significant dimen- 
sional reduction, as the aldehyde and ^A contain 65 and ~ 2200 degrees of freedom. 
Results from the six functionals are presented in Table IVII A plot of the extrapolation 
results is in Fig. El 

In general, the results are quite good. '^A[^Dexact] errors for the R groups in the training 
sets(Table IVH upper panel, off-diagonals) are small compared to both the average E^^^J^ 
(-139.60 mH) and the standard deviation in E^J^J^ (13.30 mH). Most of the extrapolations 
are also good, with '^A[^Dexact] errors that are uniformly smaller than the corresponding 
^A[^Davg] errors (diagonals of Table EH compare upper and lower panels). The ^A[^Davg] 
energy errors are fairly good, as in the variable-geometry (H-H)5 systems, providing further 
evidence that the primitive, constant-^ A functional works rather well for multiple geometries 
(see Sec HTTBl) . 



IV. DISCUSSION 



Nearsightedness and molecular similarity suggest that a rich data set of ab initio calcu- 
lations on a functional group in various small molecules contains sufficient information to 
describe the functional group's behavior in large molecules. Here we explore new methods for 
generating semiempirical electronic structure models that are parametrized to such data sets. 
In particular, we consider a semiempirical, subsystem-based model of electron correlation. 
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This model predicts the connected pair density matrix ^ A of molecular subsystems as a func- 
tional of the subsystem electron density matrix ^D. Subsystem predictions are combined 
using a previously-developed divide-and-conquer-style treatment of the atomic-orbital-basis 
^A (LRDM). The ^A functionals are used to obtain correlation-energy functionals for den- 
sity functional theory (Eq. P). The method is tested on chains of minimal-basis (H-H)5, 
which was treated as a system of four identical and overlapping (H-H)2 subsystems. The 
extrapolation abilities of ^A ^D] functionals are tested on HOC-R molecules. 

The (H-H)5 chain results demonstrate that the model works well for these simple systems. 
The results in Fig. |21 and Table lllll show that ^A \^D] functionals fitted to ah initio data 
on (H-H)2 can reproduce the (H-H)2 ^A given the correct electron density matrix ^D^^act- 
The ^A of (H-H)5 systems can be modeled quite well using four overlaid (H-H)2 ^A [^D] 
functionals, as shown by the data in Fig. El and Table IIVI Fig. |3] and Table |3 show that 
the ^A [^D] functionals work reasonably well as DFT correlation functionals. These results 
are especially encouraging given the simple, linear form of the Vcorr[^D] functionals (Eq. 
The subsystem ^A [^D] functionals can extrapolate to molecules outside of the training set, 
as demonstrated by the HOC-R results in Sec. Illl El 

An interesting finding is that dimensional reduction of subsystem ^A seems to be a 
reasonable approximation. All of the ^A [^D] functionals had significant reductions in the 
dimensionality of ^A. This suggests that real molecular environments only explore a fraction 
of the total degrees of freedom in a functional group's ^A. This dimensional reduction may 
be useful for other models of electron correlation. Our results also suggest that, for these 
systems, simple quadratic functions are a fairly good model for the input:output relation of 
the dimensionally reduced data. 

It is also interesting that subsystem ^A [^D] functionals that are defined in a basis set can 
be used for multiple subsystem geometries. For both hydrogen chains and HOC-R molecules, 
a single ^A [^D] functional provided good ^A predictions for a fairly wide range of different 
geometries. Real-space ^A [^D] functionals may be more general than those presented here. 
However, the success of the basis-set functionals is encouraging. 

The utility of our semiempirical, subsystem-based treatment of electron correlation relies 
upon its ability to provide high-accuracy treatments of systems that are too large for current 
ab initio methods. In particular, we would like to perform high-accuracy DFT calculations 
on polypeptides, using ^A [^D] or "^Dxci^D] functionals for the twenty most common amino 
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acids. In order to provide an accurate treatment of electron correlation, each of the twenty 
functional would be parametrized to high-level ab initio calculations in a large basis set 
(e.g. CCSD(T)/cc-pVTZ) 63]. These calculations would place each amino acid in a large 
number of small-molecule environments. The environments could include the amino acid in 
various 2- or 3-residue polypeptides, or with capping groups as in Ref. ^. 

As stated above, the principle challenge to reaching this goal is the computational ex- 
pense of parametrizing subsystem [^D] or "^Dxcl^D] functionals for moderately-sized 
subsystems in a large basis set. Parametrizing a [^D] functional for a group with M 
basis functions requires 0{exp{M'^)) data points. In the current proof-of-concept work, we 
studied small subsystems in a minimal basis set, so that new ^A [^D] functionals could be 
parametrized and tested in a relatively short time. Data sets for the systems studied here, 
minimal-basis (H-H)5 and HOC-R, could typically be generated in a couple of days on a 
single 2.8 GHz Xeon processor. 

Another challenge to using our approach for high-accuracy DFT on large systems is 
the computational expense of using the high-accuracy subsystem functionals, as this would 
require performing DFT calculations in a large basis set. One interesting possibility for 
mitigating this expense is to parametrize subsystem functionals that predict both and 
^A in a large basis set as a function of ^-D in a smaller basis. This would enable us to use a 
database of high-accuracy, large-basis-set subsystem calculations to correct small-basis-set 
DFT on a large system. 



We believe that the general approach presented here (see flowchart. Sec. Ill B|l may be 
useful for modeling several aspects of electronic structure in addition to intra-subsystem 
electron correlation effects. We are currently developing a semiempirical, subsystem-based 
treatment of dispersion interactions in DFT that uses semiempirical functionals to predict 
a subsystem's polarizability as a function of its ^D. Also, as mentioned above, we are 
developing a treatment of core polarization in effective core potentials. This treatment uses 
a semiempirical functional to predict the change in core electron density matrix as a function 
of the valence density matrix and the core-electron Hamiltonian. 

This work explores a new approach for taking advantage of molecular similarity in elec- 
tronic structure theory. The results suggest that it may be possible to construct accurate 
semiempirical models by extracting transferable information from ab initio data on small 
molecules. However, the applicability of the method to larger systems remains to be ex- 
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V. SUPPORTING INFORMATION 



A. Sources of error in fixed-geometry (H-H)5 

The conclusion that the error in fixed-geometry (H-H)5 is due to long-range order is 
supported by results from the (Ci = 7, C2 = 8) \^D] functional discussed in Sec IIII Bl 
This functional gave a better description of the fixed-geometry (H-H)2 systems than the 
functional in Table |n] However, this functional does not give a better description of the 
(H-H)5 systems: the D ^^act] l^^corrl is 5.01 (11.30) mH, much larger than the 1.81 (1.99) 
value in Table ITVl 

To further confirm that the fixed-geometry error is due to the effects of long-range order, 
we parametrized \^D] functionals for a new fixed-geometry (H-H)5 system with increased 
long-range order. This system was generated as in Sec. IIII Al but with an (H-H) ^ (H-H) 
spacing of 1.0 A rather than 1.6 A. Its (11-11)2 subsystems are expected to be even less 
similar to isolated (H-H) 2 molecules. The increased long-range order is seen in an increased 
(though still quite small) subsystem decomposition error, with average (standard deviation) 
2A^,„fe \6Ecorr\ of 1.01 (0.60) mH vs. the 0.02 (0.04) mH values in Table HH This system's 
^A \^D] (H-H)2 functionals gave |5i?corr| for (H-H)2 comparable to the values in Table ITTTl 
average (standard deviation) values of the Dexact\ \5Ecorr\ for the testing-set data are 
0.29 (0.43) mH. However, as expected, the increased long-range order meant that the ^ A \^D] 
functionals parametrized on isolated (H-H) 2 molecules gave very poor results when applied 
to (H-H)5. The ^A[^De2,.arf] I'^-^'corrl was 17.93 (4.89) mH, much larger than the value in 
Table HYl 

B. Training set choice calculations for (H-H)2 \^D\ functionals 

The predictions of a semiempirical model should not depend on the choice of training set 
data. We tested this by parametrizing several different (H-H)2 ^A \^D] functionals (102 for 
variable-geometry systems, 74 for fixed-geometry systems), each with a different training set 
choice. Table IVTTl presents the average and standard deviation, taken across the training set 
choices, of the average \5Ecorr\ values of each data set. 

To clarify how the results in Table IVlTl were obtained, let Schoice denote the set of Nchoice 
different choices of training set tested, where Nchoice equals 102 and 74 for variable- and 
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fixed-geometry systems. The training set choices in Schoice are indexed by x. Let Sf^^^^ 
denote the set of (H-H)2 molecules in training set x, where each molecule is indexed by ix- 
Each Sfj.^^^ contains 500 of the 1000 total (H-H)2 molecules f Table ITT|). with the remainder 
in the test set. Let the absolute correlation energy error |5ii^corr| for ^^[^Dexact] of each data 
point in training set x be denoted \6Ecorr\ i'ix), and let AVE{ } and STDEV{ } denote the 
operations of calculating the average and standard deviation of a set of points. The average 
(standard deviation) values of the first entry in Table IVlTl frow 1, column 2), denoted "A" 
and "B", are obtained as 

A = AVE{ AVE{\5E,orr\ (t.), l. e S^,,,J, X E Schoice} (10) 
B = STDEV{ AVE{\SE,orr\ (^.),^. e S^,,,J,Xe S,ho^ce} 

The results in Table I VI II verify that the [^D] functional predictions do not depend 
very much on the training set choice. 

When parametrizing a model, it is useful to test models that were parametrized with an 
incorrect input :output relation in the training set. If the model is implemented correctly, 
and is modeling a real physical relationship, scrambling the data should degrade the results. 
Table IVIll includes results from a [^D] functional where the from each molecule in 
the training set is paired randomly with the ^A of a different molecule. This scrambles the 
input:output relation of the training data, and is denoted ^A[^Dexact]{scr). As expected, this 
functional is no better (and sometimes worse) than ^A[^Davg], which uses a single choice of 
^A for all data points. 
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Approximation "^^xsub "^^PCA ^Af^De^act] 



Subsystem decomposition 


Yes 


Yes 


Yes 


Dimensional reduction 


No 


Yes 


Yes 


Prediction from 


No 


No 


Yes 



TABLE I: Types of connected pair density matrix obtained in the results, and the approxima- 
tions associated with each. 
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Variable geometry Fixed geometry 

(H-H)2 (H-H)5 (H-H)2 (H-H)5 

Ndat 1000 93 1000 99 

Point charges 4 10 4 10 

(H-H) bonds 0.5 ^ 1.0 A 0.7 A 

(H-H) ^ (H-H) 0.9 ^ 3.0 A 1.6 A 

Ci 6 4 

C2 5 5 

TABLE H: Details of the four data sets for linear dimerized hydrogen chains. N^at is the total 
number of molecules in the data set. Ci and C2 are the number of and principal components 
used in the \^D\ functionals (Eq. [T)). 
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Prediction V train V test F train F test 
"^/^[^Dexact] 1-95 (1.61) 1.91 (1.85) 0.37 (0.51) 0.46 (0.96) 
^ApcA 1.29 (1.29) 1.27 (1.22) 0.11 (0.16) 0.13 (0.27) 
2A[i£>a„e] 3.78 (3.77) 3.66 (3.47) 1.53 (2.72) 1.49 (2.35) 
MP2 85.83 (18.93) 68.52 (2.04) 

TABLE III: Absolute Ecorr error \5Ecorr \ (mH) for training and testing subsets of tlie variable- and 
fixed-geometry (H-H)2 subsystems (V and F, respectively). Values are average (standard deviation) 
across the entire training or testing set, for a single choice of training set. MP2 \5Ecorr\ values 
are included for comparison. The average and standard deviation of the correct Ecm-r values are 
-114.39 (23.43) mH for the variable-geometry (H-H)2 and -93.10 (2.56) mH for the fixed-geometry 
(H-H)2. 
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Prediction V F 

^A[^De^act] 3.43 (2.96) 1.81 (1.99) 

"^A^sub 0.02 (0.04) 0.02 (0.04) 

^ApcA 3.33 (3.16) 1.40 (1.21) 

"^Al^Davg] 9.21 (6.53) 2.55 (5.06) 

MP2 218.14 (30.06) 169.61 (4.33) 

TABLE IV: Absolute Ecorr error IdE'corrI (niH) for variable- and fixed-geometry (H-H)5 (V and 
F, respectively). Values are average (standard deviation) across the entire data set for the choice 
of training set used in Table ITTTI MP2 j^E'corrl values are included for comparison. 
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System Prediction V F 

(H-H)2 ^A[^Dhf] 5.72 (3.91) 3.21 (1.02) 
'^A[^Ddft] 4.10 (3.34) 1.24 (0.66) 
(H-H)5 ^A[^Dhf] 15.73 (5.46) 9.24 (1.48) 
^A[^Ddft] 12.02 (4.73) 4.36 (1.90) 
TABLE V: Absolute Ecorr error jJE^corrl (niH) for DFT and corrected Hartree-Fock calculations 
using [^D] correlation energy functionals {^A[^Dhf] and "^AI^Ddft], respectively). Results are 
presented for (H-H)2 and (H-H)5, variable (V) and fixed (F) geometry hydrogen chains, average 
(standard deviation) over the entire data set for a single training set choice. 
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T71 1 1 J 

Excluded 


TT 

H 


T^ 

1 


OH 


CH3 


CI 


OCH3 


H 


2.56 


1.05 


1.11 


1.17 


1.09 


1.10 


F 


1.44 


3.24 


1.19 


1.42 


1.31 


1.23 


OH 


1.33 


1.42 


1.42 


1.50 


1.21 


1.47 


CH3 


1.43 


1.27 


1.23 


2.35 


1.32 


1.30 


CI 


1.55 


1.16 


1.18 


1.57 


8.10 


1.34 


OCH3 


1.37 


1.26 


1.23 


1.47 


1.19 


1.76 


Excluded 


TT 

H 


T^ 

1 


OH 


CH3 


CI 


OCH3 


TT 

H 


6.00 


5.50 


A AO 

4.48 


A AO 

4.48 


6.62 


5.37 


F 


6.07 


5.44 


5.09 


4.80 


6.47 


4.88 


OH 


5.50 


4.36 


4.70 


4.11 


8.13 


5.05 


CHs 


6.37 


5.83 


4.76 


4.69 


6.88 


5.52 


CI 


5.01 


4.82 


4.85 


4.28 


9.09 


5.09 


OCH3 


5.67 


5.17 


4.74 


4.35 


7.70 


5.12 



TABLE VI: Absolute E^^^J^ errors (mH) for the six different HOC-R [^D] functionals. The 
rows are the results for each of the six functionals, where the R group that was excluded from 
each functional's training data is listed in the first column. The columns show the mean absolute 
^arrr" error for each of the six kinds of HOC-R molecules in the testing set. Extrapolations to the 
R group excluded from each functional are shown in boldface. Upper and lower panels are data 
for A\^Dexact] and '^/^\^Davg\- The extrapolation results for ^A[^Dea'act] are plotted in Fig. [51 
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J-^ exact i 


^Af^D 1 

■^l ^avg\ 




V train 


i. / l^U.iO ) 


0. 1 U l^U. iZ J 


D.yo (^u.ooj 


V I/cbl 


1.00 ^^U.iOJ 






V (H-H)5 


3.12 (0.23) 


10.18 (0.69) 


12.91 (0.74) 


F train 


0.45 (0.04) 


1.38 (0.11) 


1.43 (0.12) 


F test 


0.48 (0.03) 


1.36 (0.15) 


1.41 (0.08) 


F (H-H)5 


1.98 (0.17) 


2.53 (0.04) 


3.95 (0.30) 



TABLE VII: Absolute Ecorr errors |(5Ecorr-| (mH) for (H-H)2 \^D\ functional, for multiple 
choices of training set. Each entry is the average value of all molecules in the training or testing 
data set, average (standard deviation) over the training set choices (see text for details). Results 
are reported for variable- and fixed-geometry systems (respectively V and F), for (H-H)2 training 
and testing sets and extrapolation to (H-H)5 (respectively train, test, and (H-H)5). 
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Figure captions 

FIG. 1: Schematic of a nearsightedness-based divide-and-conquer treatment of electronic struc- 
ture for a generic four-element chain. The electronic structure of the three subsystems (boxed 
regions) are obtained separately (calculations I-III) and combined into an approximate electronic 

structure for the entire system ("Total"). The calculated electronic structure near the edges of 
each subsystem (dotted lines) is incorrect due to short-range edge effects, and is not used in the 
final approximate structure. 



FIG. 2: Predicted vs. real Ecorr (mH) for variable- and fixed-geometry (H-H)2 (A and B). The 
correlation coefficient between real and predicted Ecorr are in parentheses. To reduce congestion, 
the variable-geometry '^/S.pcA and /^\^Dave\ Ecorr arc shifted down by 30 and 60 mH. MP2 E^orr 
are shifted down by 86 and 67 mH for the variable- and fixed-geometry results, respectively. 

FIG. 3: Predicted vs. real correlation energies (mH) for variable- and fixed-geometry (H-H)5 (A 
and B), predicted using (H-H)2 functional. between real and predicted Ecorr are in parentheses. 
To reduce congestion, the variable-geometry '^/S.r^sub-, '^^PCA and ^A[^Da„e] Ecorr are shifted down 
by 60, 60, and 120 mH, respectively. MP2 Ecorr are shifted down by 216 and 163 mH for the 
variable- and fixed-geometry results, respectively. 

FIG. 4: Predicted vs. real correlation energies (mH) for DFT and corrected Hartrcc-Fock calcu- 
lations using (H-H)2 \^D\ functional. Results are presented for variable- and fixed-geometry 
(H-H)5 (A and B). R^ between real and predicted Ecorr are in parentheses. 
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FIG. 5: Extrapolation results. Predicted vs. real '^A[^D(,xact] E^o?r' kinds of HOC-R 

molecules. Each of the HOC-R data sets is modeled using the [^D] functional that was not 
trained on data from that R group. The correlation coefficients R^ between real and predicted 
^m?r^ are in parentheses. Absolute E^^J^^ errors for the plotted data are the diagonal (boldface) 
entries in the upper panel of Table IVII 



32 



c ) C 



( ) ( ) I 



( );( 


) ( 


); ( ) 


II 








( ) ( ) :( 


) ( ) 


III 














( ) ( 


) ( 


) ( ) 


Total 



Figure [H 

Janesko et. al. 

Journal of Chemical Physics 

33 




Figure |21 

Janesko et. al. 

Journal of Chemical Physics 



34 



I 

E 



I 

E 



LU 

■o 





-100 
-200 
-300 
-400 
-500 

-195 
-200 
-205 
-210 
-215 
-220 
-225 
-230 
-235 



'A[b3,a„] (0.99) 


□ 


(A) 


^Apcft-60 mH (0.99) 
'A,3„b-60 mH (0.99) 


+ 
o 


^Ar^D„,„]-120 mH (0.97) 






MP2-216 mH (0.98) 






1:1 - 






XI 1 1 1 







-400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 
Real E,.rr (mH) 



2a[^D,,,,J(0.83) 
^ApcA (0.99) 


□ 
-1- 




□ 


(B) _ 


2A[^D,,g](0.00) 
MP2-163mH (0.99) 


o 




o 


□ 


1:1 - 

□ 




^ 


















□ 










xsoc 




X 





-235 



-230 



-225 
Real E, 



-220 
corr »mH) 



-215 



-210 



Figure 01 

Janesko et. al. 

Journal of Chemical Physics 

35 




Figure HI 

Janesko et. al. 

Journal of Chemical Physics 



36 




-190 I 1 1 1 1 1 1 1 1 1— 

-180 -170 -160 -150 -140 -130 -120 -110 -100 



Real Bcon ("iH) 



Figure |31 

Janesko et. al. 

Journal of Chemical Physics 



37 



